#############
#MAP EMPIRES#
#############
#This to obtain:
#Figure_1a.jpg
#Figure_1b.jpg

rm(list = ls())
library(ggplot2)
library(dplyr)
library(ggrepel)
library(RColorBrewer)
library(lemon)
library(maptools)
library(ggsn)
library(ggplot2)
library(dplyr)
library(lemon)
library(raster)
library(sf)

setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")

#Reading polygons
Ottoman_1797 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="Ottoman_1797_GCS", quiet = T)

#Simplifying polygons
Ottoman_1797<-st_simplify(Ottoman_1797,  dTolerance = 0.005)


Habsburg_1797 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="Habsburg_1797_c_GCS")
#Simplifying polygons
Habsburg_1797<-st_simplify(Habsburg_1797,  dTolerance = 0.005)


square_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                   layer="polygon_GCS")
#Simplifying polygons
square_polygon<-st_simplify(square_polygon,  dTolerance = 0.005)


mil_col_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                         layer="krajna_polygon_d_GCS")
sf::sf_use_s2(FALSE)
mil_col_polygon<-st_simplify(mil_col_polygon,  dTolerance = 0.005)

study_area <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                          layer="study_area_b")
study_area<-st_simplify(study_area,  dTolerance = 0.005)

all_coun <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                          layer="all_countries_GCS")
all_coun<-st_simplify(all_coun, dTolerance = 0.0005)


euro_raster <- raster("./Dropbox/Legacies_Project/Analysis/data/euro_raster.tif")

max_lat<-52.669853 + 5
min_lat<-34.815009 - 5

max_lon<-40.447201 + 5
min_lon<-7.556495 - 5


e <- extent(min_lon, max_lon, min_lat, max_lat)
rc <- crop(euro_raster, e)	


#Lower resolution raster
rc.aggregate <- aggregate(rc, fact=5)
res(rc.aggregate)

#Get rid of really high mountains
rc <- reclassify(rc.aggregate, c(10000,Inf,NA))
rm(rc.aggregate)
rna <- reclassify(rc, cbind(NA, 0))

rm(rc)
#Convert rasters TO dataframes for plotting with ggplot
hdf <- rasterToPoints(rna); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Hill")

cols.fill = c("Ottoman Empire &\nVassal States" = alpha("blue", 0.5), 
              "Habsburg \nCivilian Territory" = alpha("pink", 0.5), 
              "Habsburg \nMilitary Territory" = alpha("red", 0.5))


map<-ggplot()+
  #Raster works but takes a long time
  geom_raster(data=hdf,aes(X,Y,alpha=Hill)) +
  scale_alpha(name = "Altitude", guide = "none")  + 
  geom_sf(data = all_coun, 
          fill=NA, color = "black", linewidth = 0.4)+
  geom_sf(data = Ottoman_1797, 
          aes(fill = "Ottoman Empire &\nVassal States"), color = "grey20", linewidth = 0.6)+
  geom_sf(data = Habsburg_1797,  
          aes(fill = "Habsburg \nCivilian Territory"), color = "grey20", linewidth = 0.6)+
  geom_sf(data = mil_col_polygon, 
          aes(fill = "Habsburg \nMilitary Territory"), color = "grey20", linewidth = 0.4)+
  geom_sf(data = square_polygon, 
          fill=NA, color = "red", linewidth = 1.5)+
  scale_fill_manual(name = "Territory", values = cols.fill) +
  coord_sf(ylim=c(34.815009,52.669853), xlim=c(9, 35))+
  theme(legend.position="left")+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  ggtitle("Habsburg and Ottoman Empires in 1797")+
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
        ggsn::scalebar(x.min = 32.5, x.max = 34.5,
           #Above you mention how long the scalebar should be
           y.min = 50.5, y.max = 52,
           #Above you mention how thick the scalebar should be
           dist = 250, dist_unit = "km",
           transform = T, model = "WGS84",
           location = "topright",
           st.size = 4,
           #Above you have the font size of the numbers below the scalebar
           st.dist =0.2,
           #Above you have the distance between the bar and the text, as a proportion of the y axis.
           height=0.2)
# #Above you have a number between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis

map<-reposition_legend(map, 'bottom left')
ggsave(map,file="./Dropbox/Legacies_Project/Paper/figures/Figure_1a.jpg", 
       height=19.42, width=20.2, units = "cm", dpi=300)


#########
#Croatia#
#########

#Reading rasters
euro_raster <- raster("./Dropbox/Education project/Data/Elevation/euro_raster.tif")
  
max_lat<-46.669853 + 3
min_lat<-42 - 3

max_lon<-19.447201 + 3
min_lon<-13.556495 - 3

e <- extent(min_lon, max_lon, min_lat, max_lat)
rc <- crop(euro_raster, e)

#Lower resolution raster
rc.aggregate <- aggregate(rc, fact=3)
res(rc.aggregate)
#plot(rc.aggregate)

#Get rid of really high mountains
rc <- reclassify(rc.aggregate, c(10000,Inf,NA))
rm(rc.aggregate)
rna <- reclassify(rc, cbind(NA, 0))

rm(rc)
#Convert rasters TO dataframes for plotting with ggplot
hdf <- rasterToPoints(rna); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Hill")


cols.fill = c("Ottoman Empire" = alpha("blue", 0.5), 
              "Habsburg \nCivilian" = alpha("pink", 0.5), 
              "Habsburg \nMilitary" = alpha("red", 0.5),
              "Study Area" = alpha("firebrick4", 0.5))

map_croatia<-ggplot()+
  #Raster works but takes a while
  geom_raster(data=hdf,aes(X,Y,alpha=Hill)) +
  #guides(fill = guide_colorbar()) 
  # use the "alpha hack"
  scale_alpha(name = "Altitude", guide = "none")  + 
  geom_sf(data = all_coun, 
          fill=NA, color = "black", linewidth = 0.4)+
  geom_sf(data = Ottoman_1797, 
          aes(fill = "Ottoman Empire"), color = "grey20", linewidth = 0.6)+
  geom_sf(data = Habsburg_1797,  
          aes(fill = "Habsburg \nCivilian"), color = "grey20", linewidth = 0.6)+
  geom_sf(data = mil_col_polygon, 
          aes(fill = "Habsburg \nMilitary"), color = "grey20", linewidth = 0.4)+
  geom_sf(data = study_area,
          aes(fill = "Study Area"), color = "grey60", linewidth = 0.7) +
  scale_fill_manual(name = "Territory", values = cols.fill) +
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201), expand = FALSE)+
  theme(legend.position="left")+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  ggtitle("Habsburg Military Area and Modern-Day Croatia")+
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
        ggsn::scalebar(x.min = 18.5, x.max = 19.1,
           #Above you mention how long the scalebar should be
           y.min = 46.2, y.max = 46.5,
           #Above you mention how thick the scalebar should be
           dist = 50, dist_unit = "km",
           transform = T, model = "WGS84",
           location = "topright",
           st.size = 4,
           #Above you have the font size of the numbers below the scalebar
           st.dist =0.2,
           #Above you have the distance between the bar and the text, as a proportion of the y axis.
           height=0.2)

map_croatia<-reposition_legend(map_croatia, 'bottom left')

#map_croatia
ggsave(map_croatia,file="./Dropbox/Legacies_Project/Paper/figures/Figure_1b.jpg", height=19.42, width=20.2, units = "cm", dpi=300)

